%
%计算有限温的格林函数
%


L=4;

[lattx, latty, latti] = ckb_latt(L);


figure (1); hold on;
for uix=1:L
    for uiy=1:L
        scatter(lattx(uix, uiy, :), latty(uix, uiy, :))
        text(lattx(uix, uiy, 1), latty(uix, uiy, 1), num2str(latti(uix, uiy, 1)))
        text(lattx(uix, uiy, 2), latty(uix, uiy, 2), num2str(latti(uix, uiy, 2)))
    end
end
%plot(1:1:total_kpt, eng2_in_path)
%plot(1:1:total_kpt, eng3_in_path)
hold off;

saveas(gca, "ckb_latt.png")


#hmlt = ckb_hmlt(L, -1, 0, 0);

#eig(hmlt)

hmlt = ckb_hmlt(L, -1, -1.0, -0.6);

eig(hmlt)

[grn, grn_t1, grn_t2] = grn_gnrl(hmlt, 1/2, 16, 0.125);
%
grn_c = -transpose(grn);
n = 0;
for idx=1:1:2*L*L
    grn_c(idx, idx) = grn_c(idx, idx) + 1;
    n = n + grn(idx, idx);
end
%计算超导关联函数
swave = 0.0;
for idx=1:1:2*L*L
    % (c^dagger_1u c^dagger_1d - c^dagger_1d c^dagger_1u) *
    % (c_id c_iu - c_iu c_id)
    % = 4 c^dagger_1u c^dagger_1d c_id ciu
    swave = swave + grn_c(1, idx) * grn_c(1, idx);
end
swave = swave% / 2 / L / L
exit

%grn
%grn_t1(:, :, 1)
%grn_t2(:, :, 1)


%超导的序参数
%第一个格子的四个最近邻
dspp = zeros(4, 3);
dspp(1, 1) = latti(1, 1, 1); dspp(1, 2) = latti(1, 1, 2); dspp(1, 3) = 1.0;
dspp(2, 1) = latti(1, 1, 1); dspp(2, 2) = latti(L, 1, 2); dspp(2, 3) = -1.0;
dspp(3, 1) = latti(1, 1, 1); dspp(3, 2) = latti(L, L, 2); dspp(3, 3) = 1.0;
dspp(4, 1) = latti(1, 1, 1); dspp(4, 2) = latti(1, L, 2); dspp(4, 3) = -1.0;

dspp2 = zeros(4, 3);
dspp2(1, 1) = latti(2, 1, 1); dspp2(1, 2) = latti(2, 1, 2); dspp2(1, 3) = 1.0;
dspp2(2, 1) = latti(2, 1, 1); dspp2(2, 2) = latti(1, 1, 2); dspp2(2, 3) = -1.0;
dspp2(3, 1) = latti(2, 1, 1); dspp2(3, 2) = latti(1, L, 2); dspp2(3, 3) = 1.0;
dspp2(4, 1) = latti(2, 1, 1); dspp2(4, 2) = latti(2, L, 2); dspp2(4, 3) = -1.0;


dspp3 = zeros(4, 3);
dspp3(1, 1) = latti(1, 1, 2); dspp3(1, 2) = latti(2, 2, 1); dspp3(1, 3) = 1.0;
dspp3(2, 1) = latti(1, 1, 2); dspp3(2, 2) = latti(1, 2, 1); dspp3(2, 3) = -1.0;
dspp3(3, 1) = latti(1, 1, 2); dspp3(3, 2) = latti(1, 1, 1); dspp3(3, 3) = 1.0;
dspp3(4, 1) = latti(1, 1, 2); dspp3(4, 2) = latti(2, 1, 1); dspp3(4, 3) = -1.0;


spp_corr(dspp, dspp, grn)
spp_corr(dspp, dspp3, grn)
spp_corr(dspp, dspp2, grn)

spp_corr(dspp, dspp, grn_t2(:, :, 40))
spp_corr(dspp, dspp3, grn_t2(:, :, 40))
spp_corr(dspp, dspp2, grn_t2(:, :, 40))


%计算配对的磁化率
sus1 = 0.0;
sus2 = 0.0;
for idx=1:1:48
    sus1 = sus1 + spp_corr(dspp, dspp, grn_t2(:, :, idx));
    sus2 = sus2 + spp_corr(dspp, dspp3, grn_t2(:, :, idx));
end
sus1 = sus1 * 6 / 48
sus2 = sus2 * 6 / 48


%计算磁关联
%以第一个元胞的格子为reference
refidx = latti(1, 1, 1);
refx = lattx(1, 1, 1);
refy = latty(1, 1, 1);
%(0,0)点和（0.5pi, 
fmsus11 = 0.0;
fmsus12 = 0.0;
afmsus11 = 0.0;
afmsus12 = 0.0;
for idx=1:1:48
for uix=1:1:L; for uiy=1:1:L
    #黑黑
    bsite = latti(uix, uiy, 1);
    bx = lattx(uix, uiy, 1);
    by = latty(uix, uiy, 1);
    #黑白
    wsite = latti(uix, uiy, 2);
    wx = lattx(uix, uiy, 2);
    wy = latty(uix, uiy, 2);
    #
    #相位
    fac = cos(0.5*pi*(refx-bx)+0.5*pi*(refy-by));
    fac2 = cos(0.5*pi*(refx-wx)+0.5*pi*(refy-wy));
    corr = grn_t2(refidx, refidx, idx)*grn_t2(bsite, bsite, idx);
    corr = corr + grn_t2(refidx, bsite, idx)*grn_t1(refidx, bsite, idx);
    corr = corr - grn_t2(refidx, refidx, idx)*grn_t2(bsite, bsite, idx);
    corr = corr - grn_t2(refidx, refidx, idx)*grn_t2(bsite, bsite, idx);
    corr = corr + grn_t2(refidx, refidx, idx)*grn_t2(bsite, bsite, idx);
    corr = corr + grn_t2(refidx, bsite, idx)*grn_t1(refidx, bsite, idx);

    fmsus11 = fmsus11 + corr;
    afmsus11 = afmsus11 + fac * corr;

end; end
end


fmsus11 = fmsus11 * 6 / 48
afmsus11 = afmsus11 * 6 / 48


%计算等时的磁关联
% c c^d
[grn, non, non2] = grn_gnrl(hmlt, 0.0, 48, 0.125);
grn_c = -transpose(grn);
n = 0;
for idx=1:1:2*L*L
    grn_c(idx, idx) = grn_c(idx, idx) + 1;
    n = n + grn(idx, idx);
end
part_number = n
afm_bb = 0.0;
fm_bb = 0.0;
afm_bw = 0.0;
fm_bw = 0.0;
%以第一个元胞的格子为reference
refidx = latti(1, 1, 1);
refx = lattx(1, 1, 1);
refy = latty(1, 1, 1);
%(0,0)点和（0.5pi, 0,5pi)
for uix=1:1:L; for uiy=1:1:L
    #黑黑
    bsite = latti(uix, uiy, 1);
    bx = lattx(uix, uiy, 1);
    by = latty(uix, uiy, 1);
    #黑白
    wsite = latti(uix, uiy, 2);
    wx = lattx(uix, uiy, 2);
    wy = latty(uix, uiy, 2);
    #相位
    fac = cos(0.5*pi*(refx-bx)+0.5*pi*(refy-by));
    fac2 = cos(0.5*pi*(refx-wx)+0.5*pi*(refy-wy));
    %
    corr = grn(refidx, refidx)*grn(bsite, bsite);
    corr = corr + grn(refidx, bsite)*grn_c(refidx, bsite);
    corr = corr - grn(refidx, refidx)*grn(bsite, bsite);
    corr = corr - grn(refidx, refidx)*grn(bsite, bsite);
    corr = corr + grn(refidx, refidx)*grn(bsite, bsite);
    corr = corr + grn(refidx, bsite)*grn_c(refidx, bsite);
    [bsite, fac, corr];
    %黑黑
    afm_bb = afm_bb + fac * corr;
    fm_bb = fm_bb + corr;
    %
    corr = grn(refidx, refidx)*grn(wsite, wsite);
    corr = corr + grn(refidx, wsite)*grn_c(refidx, wsite);
    corr = corr - grn(refidx, refidx)*grn(wsite, wsite);
    corr = corr - grn(refidx, refidx)*grn(wsite, wsite);
    corr = corr + grn(refidx, refidx)*grn(wsite, wsite);
    corr = corr + grn(refidx, wsite)*grn_c(refidx, wsite);
    %黑白
    afm_bw = afm_bw + fac2 * corr;
    fm_bw = fm_bw + corr;
end; end

afm_bb
fm_bb
afm_bw
fm_bw


